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We perform a variational quantum Monte Carlo simulation of the transition from a Bardeen- 
Cooper-Schrieffer superfluid (BCS) to a Bose-Einstein condensate (BEC) at zero temperature. The 
model Hamiltonian involves an attractive short range two body interaction and the atoms number 
2N = 330 is chosen so that, in the non-interacting limit, the ground state function corresponds 
to a closed shell configuration. The system is then characterized by the s-wave scattering length 
a of the two-particle collisions in the gas, which is varied from negative to positive values, and 
the Fermi wave number kp. Based on an extensive analysis of the s-wave two-body problem, one 
parameter variational many-body wave functions are proposed to describe the ground state of the 
interacting Fermi gas from BCS to BEC states. We exploit properties of antisymmetrized many-body 
functions to develop efficient techniques that permit variational calculations for a large number of 
particles. It is shown that a virial relation between the energy per particle and the trapping energy is 
approximately valid for —0.1 < 1/kpa < 3.4. The influence of the harmonic trap and the interaction 
potential as exhibited in two-body correlation functions is also analyzed. 

I. INTRODUCTION 

The experimental realization of a degenerate Fermi gas in 1999^, boosted theoretical and experimental efforts to 
study interacting Fermi gases, in particular, the formation of molecules and h ighly corr elated pairs from a balanced 
mixture of neutral interacting Fermi atoms in two different hyperfine spin states ^ * 4 I 5 I 6 I 7 I. The possibility of tuning the 
strength of the interaction between particles in different spin states via Feshbach resonances, results in the formation 
of Cooper pairs (molecules) for negative (positive) values of the scattering length a. At low temperatures, these pairs 
and molecules can form a Bardeen-Cooper-Schrieffer (BCS) superfluid state and a Bose-Einstein condensate (BEC) 
respectively. When crossing from the BCS to the BEC region, and viceversa, a grows in magnitude until it diverges 
at the resonance. In this limit the scattering length is no longer a relevant scale, and the properties of the gas become 
independent of the specifi c deta ils of the interaction potential. This is the so called unitarity limit in which the gas 
is assumed to be universa l 3 * 7 * 8 * 9 !, because its properties depend locally just on the density and the temperature, i.e., 
the only relevant scales in this quantum gas are the interparticle spacingand the Fermi energy. Consequently the gas 
properties can be expressed in terms of them and universal parameter d 3 * 7 * 8 ! 

Previous treatments of the BCS-BEC crossover in degenerate atomic gases have been done using different ap- 
proaches, we can mention the self -consistent many -body approach^, the effective field theory^, and more recently 
quantum Monte Carlo calculation d 11 * 12 * 13 * 14 * 15 * 16 -^. This last treatment has been predominantly based on the fixed 
node Quantum Monte Carlo technique. In most of these calculations, the two-component Fermi gas is considered as 
an homogeneous system although, experimentally, the gas has an intrinsic inhomogeneous nature provided usually by 
a magnetic and/or optical trap. Such confining can be described by a harmonic potential. An interesting parameter 
calculated in those approaches is (3, which relates the Fermi energy of the ideal Fermi gas Eifg and the total energy 
of the interacting gas E. This p arameter is expected to acquire a universal value at unitarity^. The predicted values 

for (3 ranges from -0.75 to -O^ j 7 1 10 1 11 ' 1 - 3 ^. Fi rst experimental estimates gave /3 0.3rP and (3 0.49 ± O.O^ 

while, more recently, values around —0.54^21221 have been reported. The later results are based on measurements of 
the gas cloud radii at unitarity. 

In a recent workPH, we employed variational quantum Monte Carlo techniques (VQMC) to describe a balanced 
two-component interacting gas confined in a three-dimensional harmonic potential. There, we reported direct tests 
of the universality hypothesis in the unitarity limit that include: (i) the verification of virial relations for N— 4, 10, 

20, 35, 56, 84, 120, 165 and 220, (ii) the variational estimate j3fu > —0.50^°'^' using a linear fit of the energy per 
particle. In that paper we also briefly reported an analysis on observables like the system energy and density profiles 
in the BCS-BEC crossover. In particular we found N— independent energy curve features through the crossover. 

In the present article, we extend the analysis of Ref paying special attention to exhibit additional properties of the 
trial many-body wave functions, whose structure incorporates, in an analytic and compact form, important features 
of the trapped two-body system. Particular properties of the antisymmetrized many-body functions, let us develop 
efficient techniques that permit variational calculations for an unusual large number of particles. The optimized wave 
functions allow the study of the influence of the harmonic trap and the interaction potential in energies, densities and 
two-body correlation functions, all along the crossover. The correlations between atoms in the same hyperfine state 
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show the Pauli-blocking evolution as a function of a. Similarly, the correlations for atoms in different hyperfine states 
give information on the formation of molecules and Cooper pairs. The applicability of virial relations to our results 
is also analyzed. Here we report results for N — 165 particles per each hyperfine state. 

This work is organized as follows: in section II, an extensive discussion of the two-body problem in the trap is 
done, and expressions for two-body functions that contain interaction and trap effects are obtained. The results of 
that section are then used to construct variational many-body wave functions for each region of the crossover. In 
section III, we address the many-body system and exploit the structure of the variational wave functions to optimize 
numerical calculations. There, we describe in detail the procedure for the variational quantum Monte Carlo simulation 
and energy evaluation. This section also contains the results for optimal variational parameters and energies, as well 
as densities and two-body correlation profiles. Our conclusions are presented in section IV. 

II. THE TWO-BODY PROBLEM 

In this section we shall establish the two-particle system features. As it is well known, in the limit of low energies it 
is expected that the scattering process, represented by the s-wave scattering length a, determines the general features 
of the state of two colliding particles, regardless the detailed form of the interaction potential among them. 

Here we consider two particles of mass m trapped in a harmonic potential of frequency u> and interacting through 
an isotropic attractive potential of finite range 6/2 given by 

V(r i>j ) = V e- 2 ^- r n\/\ Vo <0 (1) 

where the j and J, subindices denote two different hyperfine atomic states and b « yjh/muj. The potential is chosen 
so that, in otherwise free space, it would admit a finite number of bound states as its strength Vo is varied. 

For interactions taking place in free space, the Schrodinger equation 

+ V]<j> = £<p (2) 

has analytical s-wave solutions^ </>(r) = v(r)/r both in the continuum 

v (y) = c i J ibV£^/n(y) + ciJ^bVs^/hiv)' ( 3 ) 

and in the bound states region 

<y) = c + J b ^£^/^y) ( 4 ) 

where y = (e~ r / b , ( = (by/\Vo\m/K), and J„ represents the Bessel function of the first kind of order v. By imposing 
the proper boundary conditions and considering the limit £ — > + , the following expression is found for the s-wave 
scattering length dependent just on £ 
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with Nq the Bessel function of the second kind and order zero, and C the Euler constant. This scattering length 
diverges whenever Jo(C) = 0- Denoting the zeros of the Jo Bessel function in increasing order by zj. (k = 0, 1, 2, ...), 
the potential V(r) admits just A:-bound states for Zk < ( < Zk+i- The discrete eigenvalues are determined by the 
boundary condition at r = 0, J b ^^^/ h (b^\Vo\m/h) = 0. 

When the two-body collision process takes place in the presence of an isotropic harmonic potential, the two-body 
Schrodinger equation can be separated in a center of mass equation 

+ ~Mu 2 R 2 cm ]$(Kcm) = EcmH^cm), (6) 

and a relative coordinate equation 

[ ^ + l^r 2 + V(r)Mr) = e ( p(r). (7) 
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FIG. 1: (Color online) Lowest s-wave relative energy eigenvalues in units of htu for two colliding trapped particles, Eq. Q, 
around the first resonance. It was evaluated by considering a potential range b/2 — 0.015^/h/mu and a strength Vb starting 
from Vb ~ to the lowest |Vb| yielding a —* + . The scattering length is measured in units of yjhjrnui. 



with )jl = to/2 and M = 2m. The former is the isotropic harmonic oscillator equation whose solutions are well known, 
and the latter can be numerically solved given b and Vq. 

Figure [T] illustrates the s-wave lowest eigenvalues e„, n < 5, as a function of the inverse of the scattering length 
when the potential parameters are in the first resonance region ( £ = bW\ V^m/Ti around zq) and b « y/h/muj. 
For — oo < l/a < 0, the spectrum is discrete with positive values in counterpart to the free space system which has 
a continuum spectrum and no bound states. In fact for a — ► CP, e„ = (2n + 3/2)7kJ as expected. At resonance, 
l/a = 0, the ground state energy eo ~ \/2%uj. For positive a, eo decreases becoming zero at l/a ~ 1/2^/mw. For 
a — > + , eo takes values close to the ground state energy of the free space system, Eq. (|2), while the excited states 
energies become e„ >0 — (2n — l/2)%u>. In the whole region — oo < l/a < oo, all e„>o exhibit a similar behavior and 
are consecutively spaced among them by a factor of ~ 27kJ. In fact, at resonance, e n = (2n + l/2)7kJ. 

If C is further increased, the scattering length becomes negative until the second resonance is reached at z\. Around 
the second resonance, the eigenvalue e t +i as a function of a is similar to et around the first resonance. For instance, 
for £ = zi, the first excited state energy e\ becomes t\ = fiu/2. Meanwhile, the ground state energy eo remains 
similar to the corresponding ground state energy of Eq. ^ which decreases with growing Q. Higher values of £ yield 
analogous results, so that ek+t, t = 0,1, 2, ... as a function of l/a around the (k + l)*' i -resonance is similar to et around 
the first resonance. 

As expected, the qualitative behavior of the spectrum illustrated in Fig. 1 using the finite range interaction V, 
is in excellent agreement with the analytical results for two harmonically trapped particles interacting through a 
regularized contact potential {/IttTi 2 a / rn)8 reg {v i — rj)P2. For that problem, Busch et al found an implicit equation for 
the energy eigenvalues e, 



^ r(-e/(2M+3/4) = ^nj^ 

r(-e/(2?kj) + l/4) a ' 1 1 

and the explicit expression for the corresponding cigcnfunctions. In particular, for \a\ — > oo, the ground state energy 
is eo = hu>/2. We have checked that given a, the finite range interaction spectrum ek+t, t = 0,1,2,... around the 
(k + l)*' l -resonance reproduces with increasing accuracy the contact interaction spectrum as the potential range 
parameter b — > 0. In order to obtain such matching, shorter potential ranges b/2 are required for negative energies 
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FIG. 2: (Color online) Radial function u a {r) for interacting particles in otherwise free space. The zero-energy resonant function 
Moo(r) (solid line) tends to a nonzero constant as r — > oo, meanwhile it2.i(r) (dashed line) and uo.5g(r) (dotted line) correspond 
to increasingly bound states. Distances are measured in units of yjTi/mu). 



than for positive energies. 

From now on, we consider just short range potentials and £ around the first resonance condition. The general 
behavior of the s-ground state cigcnfunctions ^p a {f) is illustrated in Fig. 2 and Fig. 3 in terms of the functions u a 
(y a (r) = u a {r)/r) considering the free space and trapped system respectively. In both figures, the solid line represents 
Moo at resonance (a = oo). The structure of this ground state trapped wave function deviates significantly from its 
free-space analog not just at long distances both also near the origin. 

For a < 0, we have found that the numerical solution can be approximated using the following analytical compact 
representation: 

¥>«px(r) = J (zoe- r/b )e" 2 / 4?l (l + ce- 2r ' b )P{r /b) /r (9) 

where c is independent of r and P(r/b) is a polynomial function. In fact, this approximation has an accuracy higher 
than 0.01% by the proper choice of c and a fourth order polynomial P{r/b), both of which depend on Vq and b. 
The accuracy of this approximation was measured by evaluating the ratio ^p a px{f)/ l Pnum(' r ) between the analytical 
approximate expression Eq. ^ and the numerical solution. 

For £ > zq in the region of positive a, the ansatz for the ground state function is: 

Va px (r)=v(y(r))e- m " r2 / 4h g(r)/r (10) 

where v was defined in Eq. Q. The function tp apx is numerically accurate at least at the 1% level. Its structure let us 
understand the origin of the eigenvalue e ~ l/2hu>. In this case, v(y(r)) takes care of the boundary condition v(0) = 
so that the effective equation for g(r) is almost identical to its analog for the one dimensional harmonic oscillator 
without the requirement of becoming null at r = 0, thus admitting the possibility e = huj/2. 

The analytical approximations given by Eqsi 
in the study of the many-body system. 




9|T0 1 to the exact solutions of the two-body problem will be exploited 
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FIG. 3: (Color online) Radial function u a {r) for interacting particles in the presence of the trapping potential. The dot-dashed 
curve corresponds to the ground s-state for a negative scattering length U-o.e(r), the resonant function ttoo(r) is given by the 
solid curve, while it2.i(r) by the dashed one and uo.ssM by the dotted line. In this figure the wave functions have been properly 
normalized. Distances are measured in units of \/h/mu). 



III. THE MANY-BODY SYSTEM 



Let us consider the system made up of 2N fermions of mass m in two, equally populated, hyperfine states (N — 
Nf = 7V|) confined in an isotropic three-dimensional harmonic trap of frequency uj. The system is allowed to interact 
via collisions between particles of different hyperfine states. The Fermi gas is considered to be at zero temperature 
and the two-body collision process is approximated by the single-channel model described in the previous section, 
Eq. ([lj. The Hamiltonian describing such a system is: 

H = H trap + J2v{nj) 

N 2 , 2 N 

= E + (4 + r?d ] + E nk*t - *n\) (id 

(12) 



A. Variational Monte-Carlo simulations 



In a variational calculation, for a given form of the interaction potential, the optimal value of any variational 
parameter A in the wave function is determined by imposing that the expectation value of the Hamiltonian, 



Eq. (12) in our problem, to be a minimum with respect to such parameter. So that, 



For a system of 2N atoms, computing the expectation value requires the evaluation of a 6iV-dimensional integral. 
The main idea of the Monte Carlo methocP^l is not to evaluate the integrand at every one of the quadrature points, 
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but rather at only a relatively small representative sampling, where the sequence of configurations are distributed 
according to |\I/a| 2 / ( ^ a | ^ a ) ■ We use Metropolis algorithm 33 which ensures that the desired probability distribution 
is approached asymptotically. 

In this article, the generic form of the variational wave function will be different according to the region of the 
BCS-BEC crossover. Explicit details will be given for each region separately. 



1. Variational calculation for weakly interacting fermions 

First, let us consider the region in the potential parameters space where, for a given range 6/2, the amplitude of the 
potential is so small that no bound states are allowed in the homogeneous two-body problem. There, it is expected 
that the trapped ideal Fermi gas configuration gives a rough description of the system. Accordingly, a Jastrow-Slater 
wave function of the form 

*A* = ®IFG ■ F( (14) 

is assumed. Here A is a variational parameter, $ifg is the Fermi gas wave function given by the product of Slater 
determinants (one for each hyperfine state) describing a noninteracting system of harmonically trapped atoms, and 
the Jastrow function will explicitly include the effects of the interaction potential. It is expected that this kind of 
wave function gives an appropriate description of the weakly interacting Fermi gas in the normal regime but not in 
the superfluid one. 

The inputs of the Slater determinants are the single-particle eigenstates of a non-interacting particle in a harmonic 
trap, </>^°(r), with quantum numbers n at the position r. This construction ensures that the wave function is totally 
antisymmetric under the exchange of identical atoms. The energy of each single particle state is characterized by 
three integer quantum numbers n= (n x , n y , n z ): 

3 

E n = 7kj(- + n x + n y + n z ) , (n, = 0, 1, 2, . . .) . (15) 
Writing n = n x + n y + n z , the degeneracy of each energy level is (n + l)(n + 2) /2. A typical basis state has the form: 

€° n „ W = ( ^V /4 TT g^^ e-tV**., (16) 

110 £= x,y,z V s 

where 7J Tt? (^/ 'a^o) are the Hermite functions of order and a^o — y/K/mu). In this paper, we consider closed shell 
configurations so that the ground state is built up by taking all single-particle states with energies increasing from 
Eq = 3huj/2 up to the Fermi energy Ep = {M p-\-"S/2)huj, where Mp is the maximum energy level for a given number 
of particles. For large N, Ep ~ (6iV) 1 / 3 ?ia; and the corresponding radius is Rp — 2Ep/muj 2 . 

In the literature of interacting bosons and fermions, the Jastrow wave function usually takes the form of a product 
EL j fij °f correlation functions / that depend on the degrees of freedom of the pair i, j of interacting particles. In 
RefsP^2U!ini3 / is a function of the interparticle distance r that solves the free-space interacting two-body problem 
up to a healing distance d after which it is restricted to become constant. In those works, the parameter d is chosen 
by minimizing the energy. 

In this paper, we shall consider trial many-body wave functions which yield a continuous Fj? and continuous 
derivatives; the optimal variational parameter A of the trial wave function for the many-body system will establish 
an effective d as we illustrate below. In fact, we have studied two options for the Jastrow function: 

(i)/y = cxp(-A./iVoe" 2ri <J /b ), so that, 



Fin = exp[-A.n E y d r *T - ^il)] ( 17 ) 

(ii) 

fij = J (z e- r - /A -)(l + ce- ar «/ A «)P(r w /Aja)/rij (18) 



The first choice of the variational wave function (17) has the advantage of becoming exact when no interactions 
between hyperfine states are allowed (Aji — 0) which is the trapped ideal Fermi gas limit, where the only correl ations 
are those imposed by the Pauli exclusion principle. It is inspired on previous calculations for the nuclear matter^E!] 
where an appropriate choice of the potential allows to explore dynamically the interplay of the nuclear-to-quark matter 
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regime. In addition, this form of the variational wave function allows to estimate the energy expectation value by 
computing only spatial dependent functions in a Monte Carlo simulation^, with no need of calculating the spatial 
derivatives of the trial wave function, as we show below. 

The second choice is inspired on the general structure of the two-body wave functions in free space at low energies, 



Eq. (|9|. It allows to numerically explore shorter potential ranges than the first option ( 17 1. It reproduces the fact, first 
noticed in BCS theories, that even the slightest interaction can lead to two-body long-range-correlations, implicit in the 
polynomial P(r ij /6). Besides, it increases the reproducibility of interaction effects at short interparticle separations 
through the factor proportional to c. Deviations from Aj2 = & should be interpreted as a many-body effect. 

The structure of the variational wave function for the BCS region, allows to simplify the expectation value of the 
kinetic energy operator through an integration by partiP^, so that 

Wxa^7) - ElFG + 2 h ,hi *" ' 

where Eifg is the energy of the 2N non-interacting trapped Fermi atoms. For closed shell configurations, this energy 
can be computed using the following equation^ 



Eifg 3 



2N 4 



M l 



Huj, (20) 



instead of the large N limit, Eifg/^N = iAApfioj /4, which produces a slightly underestimated value. Eq. ( 20 1 is 



valid in general. The extra term in Eq. (19 1 reflects the increase in the kinetic energy of the system, relative to the 
Fermi-gas estimate, due to interactions. 



In the case of Eq. (17), we define the factors Wa 71 and V\ J1 through the equations 



N N 2 2 2 N N 

E E ^(^IViCiog/yj-ViOog^)!^) = ^7-E^il E ViVfajj-ViVfajOl^x) 

i=1 jj'=i i=1 i,i'=i 



— \2 



and 



ij 



XjiWxniVxnlVxjr), (21) 



(^\V(r^)\9 Xjl ) = VxjA^XjA^Xj,), (22) 



so that the expectation value of the total energy is: 

E(\. n ) = E IFG + 2A 2 71 W Ajl + V Ajl . (23) 

The two functions that remain to be evaluated (W\ J1 and V A/1 ) are local; their expectation values may be computed 
via Monte Carlo techniques as described above. A similar approach can be used in the case corresponding to the 



Jastrow function Eq. (18) 



We have performed calculations of the energy for a fixed value of N, the range b/2 and the scattering length a, 
exploring for several values of the variational parameter, picking up the one which minimizes the energy. Each run 
used about 10 3 steps for thermalization and about 10 4 more to take data. In the first rows of Table I, we rep ort the 
numerical optimal energies using the first choice for the Jastrow function and a potential range b/2 = 0M5y/hJ 
Similar results are obtained when the second choice of F J is used. The data corresponds to TV — 165, which fills eight 
shells (A4 f =8) for the harmonic potential in three dimensions; kp represents the Fermi wave number associated 
to the ideal Fermi energy Ep — (frkp) 2 /2m. For N = 165 the corresponding energy per particle for an ideal Fermi 
gas is Eifc/2N = 7.5huj, while kp = {2M F + 3) 1 / 2 ^Jmui/Ti ~ 4.3589 yjrrujjjh. The quoted error bars take into 
account the minimization process itself as well as effects of the initial conditions that could not be erased during the 
thermalization process. It is important to point out that the variational energy for the highest l/kpa value coincides 
with that obtained from a perturbative calculation using a contact interaction as can be verified from expressions 
obtained in Ref.^. 

The optimal value of the variational parameter determines the shape of the Jastrow cor rela tion function. As an 
illustration, in Fig. UJwe plot the behavior of the optimal two particle function of Eq. |l7| for b = 0.03 yjTi/muj 
and two different values of the scattering length l/kpa = —0.3873 and —15.8888 respectively We observe that the 
distance at which fij ~ 1 is larger than the potential range, this suggests long distance correlated pairs, as expected 
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FIG. 4: (Color online) First choice Jastrow correlation function for b = dashed and solid lines correspond to 

l/fcj?a = —0.3873 and —15.8888 respectively. Distances are measured in units of y/%Jmuj. 



for a BCS-like pairs. This result is also in good agreement with previous findings reported id 11 * 26 ! where the healing 
distance is used as a variational parameter. 



As a approaches the crossover region it is expected that the trial function Eq. ( 14 I will not describe properly the 



interatomic correlations: pairing effects become essential so that the quantum numbers in the Slater determinants in 
are not representative of the physical situation. 



2. Variational calculation on the BCS-BEC crossover region 

In a theory originally put forth by Eagles^ and later by LeggettP^, it was proposed that a BCS wave function of 
the form 

^ XEL =A[<t>(l h l l m2 h 2 l )...ct>(N h N l )}, (24) 

with A the antisymmetrizer operator that ensures the correct properties under particle exchanges, was more generally 
applicable than just to the weakly interacting limitP^ a BCS-like wave function could eventually describe the ground 
state from a Cooper pairing region to a BEC of composite bosons made up of two fermions. 

Following this point of view, here we propose a family of single-parameter variational wave functions for the BCS- 
BEC crossover regime, taking <fi(if, jy) as a variational extrapolation of the ground state solution of the trapped two 
body problem 

Hr lh v 3l ) £* tpin^e-^lm+ml 2 /*. (25 ) 



The variational parameter Xel modulates the optimal shape of the cloud. The wave function (24 1 using the basis (25 1 
guarantees that the Monte Carlo dynamics will be guided by effects of both paired-particles relative r,^ = — 
and center of mass Rjj = (r^ + ij-jJ/2 vectors. 

It is worth to mention that at difference with previous calculations^^ here, by explicitly including easy inter- 
pretable inhomogeneous features in the wave function, we are able to explore the trapped atoms as a whole as they 
evolve into the interacting regime. Besides, at difference with the mean field approach, no optimal individual particle 
wave functions are searched, but the global effect of the interaction on the paired-particles wave function. 

To estimate the energy in the BEC side, the algorithm described for the BCS region is not useful because it depends 
on the explicit structure of the wave function, written in a Jastrow-Salter form. In order to set the variational energy 
in a form suitable for Monte Carlo estimations we exploit the two-body structure of the potential and the primitive 



wave functions <j). The antisymmctrized wave function (24) can be explicitly written as 



A' 



*A M =£(- i rn < ^' :p w)> (26) 



where the summation is taken over all possible permutations V on set |, and cj)(i,V(i)) are wave functions having the 
form of Eq. (25) and argument (r^, r-p(^). We can split the Hamiltonian of the system in a pair-like sum, using the 
center of mass and relative coordinates of possible pairs as: 
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N 2 

z 

4-1/ fV U P ^o(») , M , ,2 p2 



with "Po any given permutation. 
Equation (25 1 let us write: 



(27) 



Ne + N 



3fru>\ 



EL 



a - ad ]>>ir n w> no) 



N 



i,V 



1 = 1 



N 



^(-1) P 2 



(28) 



with eo the ground state eigenvalue of the two body-problem. To evaluate the last two terms via a Monte Carlo 
simulation we proceed to complete the potential by adding and subtracting the term used in the two-body solution, 
then: 



Neo + N'-^ + ^Vin,) 



N 



i,V Z#i 



M 



2 R 2 

2 -- n i,n*) 



which can also be written in terms of the minors Ci a (^f \ EL ) associated to the 

N 



(29) 



(30) 



where <\>^ a represents any of the z-row wave functions. 

r 3hu}\ E L 



[iVe 



N- 



En- 



M 

E c i,«' 0--^EL)^ 2 Rl a -V{r i<a ) <f>(i,a) 



This expression results quite convenient for the simulations to be performed, where one can also take advantage from 
the relation between minors and the elements of the inverse of the transposed matrix^, 



<t>i,a = (<P ) ia 



(31) 



Thus, we can write 



Neo + N'-^+^V^) 



N 



^2<pi, a <pi, a (1 - A| L )^w 2 i? 2 ? Q - V(r ha ) 



(32) 
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As in the BCS calculation, we can sample the system using a Metropolis-Monte Carlo algorithm and estimate the 
energy as a function of the variational parameter. 

In Table I, we illustrate the results on the optimal variational parameter and different energy estimates for several 
scattering lengths. The strength of the potential was taken in the region around the first zero-res onance condition 
(zo/b) 2 Tiuj = %. The upper set of results were obtained using the Jastrow-Slater wave function Eqs. ( 14|l7 l. All other 



results considered a wave function of the Eagles-Leggett form, Eq. (24 1, with the two-body functions Eq. ( p5| taking 



(p{r^j) as the approximate solutions of the two-body problem for a given scattering length, Eqs.(9flT0 1. 

In the reported calculations using Eagles-Leggett wave functions, the range of the potential was taken as b/2 = 
0.00375^/V mw before the unitarity limit, l/kpa — 0, and b/2 = 0. 0025 Wh/mui for kpa > 0. Actually, calculations 
were performed for several potential ranges b/2 all over the crossover. The b/2 ranges reported in this table are the 
shortest for which reliable numerical results were obtained. 



For l/kpa < —0.45 the variational energy E/2N for the wave-function ( 24p5| is higher than that obtained with 



the Jastrow-Slater trial wave-function, while for l/kpa > —0.45 the situation is inverted and the BCS-wave function 

gives a lower upper bound for E/2N . A similar effect has been found in RefP^for the homogeneous gas. Beyond the 

unitarity region, i. e. for a > 0, the contribution cq/2 coming from the trapped ground state two-body eigenvalue has 

been subtracted. For 4 < l/kpa < 12 an optimal value of Xel ~ 1 yielding a local minimum was found. However, for 

Xel > 1 the corresponding mean value of the energy can be made arbitrarily small by considering Xpp large enough. 

This variational instability is expected at the extreme BEC region for any attractive potential of finite range as 

discussed previously in RefP^lfor a homogeneous gas. It could eventually be avoided by adding a repulsive interaction 

at distances much smaller than the range b/2 as suggested in the original work by LeggettP^. Implementing this idea 

within our numerical approach is very difficult since we have already set b <C y/h/muj. Although a local minimum 

was found for 1/fc/a > 4, finite range effects are expected to be significant on the reported data. 

The wave function having Xpp ~ 1 is expected for a molecular gas when Pauli Blocking effects between the 

constituting fermionic atoms are approximately compensated by the attractive finite range interaction effects. If the 

atoms did not interact, the corresponding energy per atom would be E/2N — l.haw corresponding to an ideal gas 

of trapped Bose molecules. It is also important to emphasize that, for a contact interaction, molecules formed by 

Fermi atoms are expected to have a weakly repulsive interaction, with a molecule-molecule scattering length given by 
a -npJSl 



3. Virial relations 



In the fifth column of Table I, we also report (mu> 2 R 2 ), that is, twice the mean value of the trapping potential 
energy per particle, which is more feasible of experimental verificatiorPSl than the total energy E. Notice that in the 
crossover region with —0.1 < l/kpa < 1.4, the energies in the third column are similar to (mu) 2 R 2 ). 

At unitarity, l/kpa = , a virial relation of the form E/2N = (muj 2 R 2 } is expected from previous experimental 



and theoretical studies^Sl. During the revision process of the present article, viri al the orems for trapped inter acting 
atoms outside the unitarity limit have been established both at finite temperatur e! 40 * 41 ! and at zero temperatur e 1 42 ! 43 ! 
The latter considered several forms of the interaction potential. In particular for a contact interaction with a strength 
determined by the scattering length a, it has been shown that: 

E , 2p2 > 1, dE/2N 

— = (mw R ) - ^k F a— (33 

2J\ 2 okpa 

= (m ^ )+ l^. (34) 
2/cfo al/kpa 

For a contact interaction in the free space, the energy of the bound state is eg = —h 2 /ma 2 . So that, kpade^j dkpa — 
— 2eg and, in the BEC side of the crossover, 

[mu R ) virial - [E/ZM - e /2) - - — — — . (35) 



2 



kpa 81/ kpa 



Although we have made all calculations with a finite range potential we would like to evaluate how compatible our 
results are with those arising from contact interaction predictions. Thus, in the last column of Table I we report the 



mean value of twice the potential energy per particle associated to the trap, evaluated numerically using Eq. (34) in 
the BCS side and the following expression in the BEC side: 

{mu 2 R 2 ) mnal = (E/2N - e b j2) - \ \^- W^f^l ] (36) 

2 Ikpa ol/kpa J 



11 



1/kpa 




E/2N 


eo 


< muj 2 R 2 > 


< mUJ 2 R 2 > v irial 




±0.05 


[Tilo] 


[hco] 


[hu] 


[Hlu] 


-15.88895 


-0.6 


t A f r\ 1 r\ r\r\ a 

7.450± 0.004 


1.49 


7.5 ±0.2 


<-? a -\ d i r\ r\r\f 

7.418± 0.006 


a cnooo 

-9.60883 


-0.6 


*7 ^OK_l_ A AAK 

A42o± O.OOo 


1 4 O 

1.48 


AO ±0.2 


^7 QA_L A AATK 
1 .3»± 0.00 (0 


a c\ f c o c\ 

-4.26632 


-0.6 


7.374±0.009 


1.46 


7.53 ±0.2 


7.3l± 0.03 


-2.02479 


-0.9 


7.345±0.012 


1.42 


7.53 ±0.2 


7.25±0.03 


1/kpa 


\°*lh/rrw 


E/2N 


eo 


< mw 2 i? 2 > 


< m,UJ 2 R 2 > v irial 




±0.005 


[hu] 


[fiw] 




[hu] 


-0.43893 


0.142 


7.07±0.06 


1.15 


7.72 a ±0.2 


6.88±0.09 


-0.22418 


0.134 


6.65±0.06 


0.99 


7.57 ±0.2 


6.42± 0.09 


-0.100 1 1 


A 1 £?A 

0.160 


6.3/±0.06 


A *7A 
0. /9 


6.o8±0.2 


a Ol _1_A AA 

6.21±0.09 


-0.03745 


0.182 


6.10±0.06 


0.72 


5.64±0.2 


5.88± 0.09 





0.186 


5.25±0.08 


0.50 


5.32±0.2 


5.25±0.12 


1/kpa 




(S/27V) - e /2 


— eo 


< mw 2 i? 2 > 


< mLU 2 R' 2 > v i r ial 




±0.005 


[Tito] 


[nw] 






0.13959 


0.190 


4.78±0.07 


0.21 


4.92±0.2 


4.85±0.11 


0.34876 


0.191 


4.18± 0.07 


2.45 


4.50±0.2 


4.52±0.11 


0.69684 


0.256 


3.57 ±0.07 


9.95 


3.75±0.2 


3.95±0.2 


1.04427 


0.270 


3.35±0.06 


22.65 


3.51±0.2 


3.85±0.2 


1.39107 


0.380 


3.14±0.08 


40.76 


2.99±0.2 


3.69±0.2 


2.08293 


0.665 


2.60±0.08 


93.96 


2.95±0.2 


3.38±0.2 


2.77266 


0.90 


2.0±0.1 


171.157 


2.89±0.2 


3.09±0.2 


3.46055 


0.94 


1.38±0.15 


274.99 


2.59±0.25 


2.45±0.23 


4.1469 


0.99 


1.1±0.15 


404.64 


2.59±0.25 


1.71±0.23 


5.51634 


0.99 


0.95±0.15 


756.643 


2.59±0.40 


1.19±0.23 



"The evaluation of mean radii for extended clouds requires special care of the statistics. 



TABLE I: Optimal variational parameter A, energy per particle, two-body ground eo energies, mean value of the trap potential 
energy per par ticle < muj 2 R 2 > from Monte Carlo calculations and < muj 2 R 2 > V i r iai evaluated using the virial relation given 
by Eq. ( |34|36| |. All of them were calculated as a function of 1/kpa considering 2N = 330 particles. The upper set of results 
used the Jastrow-Slater wave function Eqs. ( 14|17 1. All other results considered a wave function of the Eagles-Leggett form, 
Eq. p4| , with the two-body functions ¥?(rj,j) taken as the approximate solutions of the two-body problem for a given scattering 
length. 



with e[j the two particle energy for a finite range interaction state in otherwise free space. The numerical evaluation 
of the derivative was preceded by a numerical smoothing of data. We observe that, although the calculations were 
performed using a finite range potential, there is a reasonable agreement between the trap energies and the virial 
expressions Eqs. (34) and (36l for —0.1 < 1/kpa < 3.5. 

It was also found that as 1/kpa — > CP and 1/kpa — > + the derivatives in Eq. ( 34 1 and Eq. ( 36 > respectively attain 
a minimum. This minimum together with the small difference between the value of eo/2 and e^/2 compared to E/2N 
for < 1/kpa < 1.5, let us understand the observed similarities between the third, fifth and sixth columns of Table I 
for -0.1 < 1/kpa < 1.4. 

On the BCS side of the crossover, for —0.5 < 1/kpa < —0.2, our wave functions yield a 15% higher trapping energy 
than the virial relation predicts. It is important to mention that, for these values of 1/kpa, the atomic cloud is quite 
extended and the evaluation of both the mean energy E/2N and the mean square radius (R 2 ) requires special care 
of the statistics sampling. Improving the form of the variational wave function in this region, could diminish the 
discrepancy with the virial relation. Notice that, in the language of BCS theory, these region delimits the transition 
of the atomic cloud from a normal to a superfluid state. 

In the region 1/kpa > 4 there is also a discrepancy with the contact virial relation; it could be due to finite range 
potential effects as expected from the variational instability reported above. In fact, for these scattering lengths the 
difference between atoms interacting through a contact and a finite range potential is already evident by comparing 
their corresponding two-body ground state energies eo- These energies are in general similar but, as expected, the 
bigger differences appear for the deeply bound two-body states, a — » + . For instance, at the bottom of Table I, 
eo = — 756.6?iijj in contrast to the solution of Eq. ^ which yields -578.1S.W; for the other states the difference is less 
than 10% up to 1/kpa < 2.5 and around 15% for the remaining reported data. 
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FIG. 5: (Color online) s-ground state variational energy per atom E/2N and {muP'R 2 ) for trapped particles at unitarity as a 
function of the Fermi number Mf for closed shells. Energy is measured in units of Twj. 



4- Unitarity 



At unitarity, \a\ — ► 00, we have estimated the energy using a variational wave function of the form Eq. (26 1. 
The numerical results were presented and broadly discussed in RefJ^S in connection with the universality hypothesis. 
Notice that, in Table I, the numerical errors at unitarity are slightly larger than those obtained for nearby scattering 
lengths. The reason can be traced back to the qualitative difference between the two-body wave function Uij = ripij 
determined by Eq. (10 1. The high derealization of the unitarity paired-atoms wave function makes more difficult 
the evaluation of the energy expectation value at this limit. From these results, an upper bound to the universal 
parameter 8, defined by Ejj = EjfgV^- + ~0> is found to be p = —0.51 ± 0.01. 

In Fig. p] we show the mean value < muj 2 B 2 {M.f) > together with E{M.f) for closed shells with A4jr < 9 
corresponding to N = 4, 10, 20,35, 56, 84, 120, 165 and 220 particles. No significant discrepancy among these mean 
values is observed and the virial relation is thus verified. As reported in Ref.^, a linear relationship between the 
energy per particle at unitarity and the shell number M.jr is also found: 



E V /2N ~ (0.53 ± 0.01)(7W F + (1.95 ± 0.06))^, 



(37) 



when this expression is compared with the ideal Fermi gas result, Eq. (20 1, one obtains an upper bound for the 
universal parameter p = — 0.50^'^. 



B. 



Densities and correlations 



The information encoded in the single-particle and the two-particle correlation functions is important since those 
functions reflect the quantum mechanical nature of the particles and their collective behavior, driven by the interaction 
and trapping potentials. 

In the following, we illustrate these correlation functions for N = 165 using the optimized wave functions in the 
BCS (a < 0), the unitarity (a — > 00) and the molecular (a > 0) regimes. We have chosen examples in the crossover 
with |l/fci?a| < 1 due to its expected independence on the details of the calculation. 
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FIG. 6: (Color online) Normalized correlation function g(r) for particles in the same hyperfine state as a function of their 
relative distance. Square, circle, and triangle symbols correspond to the ideal Fermi gas, BCS (l/fc^a = —0.224) and BEC 
(1/kpa — 0.697) respectively. Distances are measured in units of yfhjmuj. 



For completeness, let us recall that the single-particle correlation function, that is, the density profile as a function 
of the distance to the center of the harmonic trap, has already been illustrated in Fig. 2 of Ref. 2 -. There, we saw 
that the trap effect is reflected by decreasing the particle density until vanishing around the Fermi radius, i.e., the 
inhomogeneous environment created by the harmonic confinement affects all the regimes as it is already evident for 
an ideal Fermi gas in the Thomas- Fermi approximation^. The shape of the BCS density profile is similar to the one 
corresponding to the ideal gas but with a different mean radius. The density increases at the center while decreases as 
it goes to the edge of the trap. These deviations can be attributed to the optimal value of the variational parameter 
which captures the interaction and correlation effects in the many-body system. This kind of shape prevails up to 
the unitarity limit. The major differences in the particle density for each regime occur around the center of the trap, 
particularly for the BEC regime where most of the paired atoms are located near the origin. 

In order to exhibit the quantum behavior of the fermionic atoms, the two-particle correlation function for particles 
in the same hyperfine state, g(r\ was computed. The calculations involved finding the fraction of atoms in the same 
hyperfine state within a relative distance (r, r + dr), as generated by the Monte Carlo sampling, irrespective of the 
center of mass position; g(r) was normalized dividing by N(N — l)/2 to account for the combinatorial of the atoms. 
Figure [6] illustrates the resulting correlation functions. The Jastrow-Slater wave function in the limit of an ideal 
Fermi gas (Xji — 0), exhibits the Pauli blocking arising from the fermionic nature of the atoms. The BCS trial wave 
function shows an slightly diminished Pauli blocking for short distances. In the molecular side, it is observed that 
particles in the same hyperfine state can be found around the same region. Although at the deep BEC regime Pauli 
blocking still inhibits the presence of atoms in the same hyperfine state, the radius at which it is evident becomes 
very short. As a consequence, if an exclusively attractive interacting potential is considered and it is large enough, 
Pauli blocking is not able to avoid a variational collapse as discussed above. All of these correlations decrease for long 
relative distances as a consequence of the presence of the trap. 

The two-particle correlation function for atoms in different hyperfine states was computed in RefPS. There (Fig. 3) 
its behavior in the BEC regime is compared with respect to the ideal regime, as a function of the relative distance 
fi,j = |r*t — r jl\ among them. It was evaluated in a similar way to g(r), taking care of the proper normalization 
factor (^V 2 ) and keeping the information of the center of mass position of the pairs. Molecule formation was indicated 
by the increase in the correlation for very short distances, r^j <C y/h/mu. Most molecules are formed for R cm < 
1.09 ^/V 

mu. An enhancement of the probability of finding pairs of particles separated at relative distances of the 
order of r^j ~ y/fi/muj indicated molecular condensation effects. 

Here, in Fig. [7] we illustrate the differences between the two-particle correlation functions of atoms in different 
hyperfine states, AKfrij, R cm ) for the ideal and BCS regimes. As in Fig. 3 of RefpH, it shows results for a set of 
radius R cm measured from the center of the trap. We observe that, although not zero, the difference is very small 
(see the abscissa scale) compared to the result for the BEC regime, in addition strong oscillations in r^- are seen for 

all Rem- 
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FIG. 7: (Color online) Probability difference AK(rij, R cm ) that two particles with different spin are found separated a distance 
nj in the BCS and ideal regimes. Each curve in this figure correspond to a spherical radius i? cm measured from the center of 
the trap. Calculations are performed at 1/kFd = —0.224. Distances are measured in units of y^h/mco. 



IV. CONCLUSIONS 



We have studied an interacting two-component Fermi gas confined in an isotropic harmonic potential in three 
dimensions. To be specific, we investigated the transition from a Bardeen-Cooper-Schrieffer state to a Bose-Einstein 
condensate at zero temperature for a system composed of N — 165 particles of equal mass in each spin-state. The 
interaction between particles of different spin was considered to be an attractive potential with very short range 
interaction; under such conditions it is expected that the many-body ground state depends just on the product of the 
scattering length a and the Fermi wave number fcj?. The BCS-BEC transition was followed as a function of kpa. 

To model the gas, we proposed a family of many-body trial wave functions for the BCS (a < 0) and the BEC 
(a > 0) sides. For small negative values of a we described the atomic gas by a Jastrow-Slater wave function. While, 
for other values of a, following Eagles and Lcggett proposal, a wave function written as the antisymmetric product of 
two-particle states was used. The two-body basis was formed by analytical compact functions that contain collision 
and trapping effects. For a given interaction range, using variational Quantum Monte Carlo simulations, we found 
the variational parameter X opt that minimizes the energy per particle of the whole system. Efficient algorithms 
to estimate the energy expectation value, exploiting properties of the antisymmetrized many-body functions, were 
elaborated to perform calculations for such a large number of particles. After considering several values of the range 
of the potential and studying the stability of the results, we reported the numerical data corresponding to the lowest 
value of the potential range that gave reliable numerical results. The corresponding optimal variational wave functions 
lead to predictions for the main properties of the trapped system like energies, mean squared radii, one and two-point 
correlation functions. 

The system energy was computed all along the crossover, and at unitarity an upper bound to the universal parameter 
j3 is found to be (3 = — 0.51±0.01. This result is compatible with the result reported in Ref.^where (3fu = — 0.50^q'q4 
was found by comparing Ejj and Eifg for M.p < 9. Those calculations indicate that the universal hypothesis yields 
results consistent with theoretical calculations even for a small N. So that, for zero temperature, the energy of a 
balanced mixture of interacting trapped fermions has the form Ejj ~ l/2(Aij? + 2)2NTlu> similar to the ideal Fermi 
gas equation Ejpq — 3/4(Aijr + 2)2N~hoj. In addition it was shown that not only at unitarity but also over the 
crossover region —0.1 < l/kpa < 1.4 the mean value of the atomic gas squared radius can be used to give a rough 
estimate of the energy per particle, since < moj 2 R 2 >^ E/2N — £q/2, where eo is the two-body ground state energy 
eo for trapped fermions for a > and zero for a < and \a\ — * oo . 

In agreement with previous results^, the energy function E(l/kpa) along the crossover follows a curve that properly 
normalized is independent of M.p. By evaluating its numerical derivative a minimum was found at unitarity . Besides, 
this function was shown to satisfy a virial relation in the interval —0.1 < l/kpa < 3.4. This relation was based on 
virial theorems for trapped atoms developed by other groups in the last few monthj 4 - * 41 * 42 * 43 !. 

The starting point of the crossover from the BCS side could be regarded as the value of l/kpa for which the 
antisymmetric product of two-particle states gives a lower expectation value of the energy with respect to the Jastrow- 
Slater wave function. According to our calculations this already occurs at l/kpa ~ —0.45. This value is similar to 
that at which a variational calculation based on scaled antisymmetric product of harmonic oscillator wave functions 
can not be applied since no minima exists 3 -^. Although the Eagles-Leggett wave function, built using the solutions of 
the two-body problem, gives a better description of the system than its Jastrow-Slater analog, it yields a mean radius 
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for the atomic cloud 15% larger than that predicted by the virial relations. Thus, it would be important to work out 
an improved trial wave function to describe this transition zone. 

In the extreme BEC region there is a variational instability which arises from the usage of finite range attractive 
potentials between the fermions. In our calculations the extreme BEC region starts when the variational parameter 
Xel yields a local minimum energy for A^J > 1- In such a case the effect of Pauli blocking is supersede by the very 
strong short range attractive potential. For 12 > l/kpa > 3.5 and 6/2 = 0. 0025 yjfijrnuj a local minimum was found 
with Xel ~ 1- This value of Xel corresponds to a many-body wave function of an ideal Bose gas of trapp ed mol ecules. 
This function does not satisfy the virial relation for a contact interaction so that even though 6/2 << ^/hj moj, finite 
range effects are not negligible for those values of kpa. 

Finally, we calculated the one-particle and the two-particle correlation functions for the BCS and BEC regimes 
and for the unitary limit. The results show that the correlation length between pairs can be much larger than the 
interaction potential range. As expected, the inhomogencous environment resulting from the harmonic confinement 
affects all the regimes. We obs erve that in the BCS regime, the paired atoms have a large correlation length particularly 
for 0.6 < R cm < 1.2-y/Ti/ mux. Pauli blocking effects were also sensible to trapping and interaction strength. Thus, 
we conclude that the approximate analytical wave function used to describe the trapped interacting gas gives a good 
compact representation of the system through the crossover region. 
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